sapply(c(
         "rjson", 
         "data.table", 
         "dplyr", 
         "ggplot2", 
         "stringr", 
         "purrr", 
         "foreach", 
         "doParallel", 
         "patchwork", 
         "lme4", 
         "lmerTest",
         "testit",
         "ggpubr",
         "latex2exp"
         ), 
       require, character=TRUE)
## Loading required package: rjson
## Loading required package: data.table
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: ggplot2
## Loading required package: stringr
## Loading required package: purrr
## 
## Attaching package: 'purrr'
## The following object is masked from 'package:data.table':
## 
##     transpose
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## Loading required package: doParallel
## Loading required package: iterators
## Loading required package: parallel
## Loading required package: patchwork
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: lmerTest
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
## Loading required package: testit
## Loading required package: ggpubr
## Loading required package: latex2exp
##      rjson data.table      dplyr    ggplot2    stringr      purrr    foreach 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE 
## doParallel  patchwork       lme4   lmerTest     testit     ggpubr  latex2exp 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE
sf <- function() sapply(paste0("./Functions/", list.files("./Functions/", recursive=TRUE)), source) # Source all fxs
sf()
##         ./Functions/General/FindDelay.R ./Functions/General/FindDelayAlt.R
## value   ?                               ?                                 
## visible FALSE                           FALSE                             
##         ./Functions/General/Utils.R ./Functions/Model/Models/RunMBayesLearner.R
## value   ?                           ?                                          
## visible FALSE                       FALSE                                      
##         ./Functions/Model/Models/RunMBayesLearnerDiffBeta.R
## value   ?                                                  
## visible FALSE                                              
##         ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayes.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayesAlt.R
## value   ?                                                                
## visible FALSE                                                            
##         ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayes.R
## value   ?                                                          
## visible FALSE                                                      
##         ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayesAlt.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/Models/RunMQLearner.R
## value   ?                                      
## visible FALSE                                  
##         ./Functions/Model/Models/RunMQLearnerDecayTo0Inits.R
## value   ?                                                   
## visible FALSE                                               
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPrior.R
## value   ?                                                      
## visible FALSE                                                  
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffBeta.R
## value   ?                                                              
## visible FALSE                                                          
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffLR.R
## value   ?                                                            
## visible FALSE                                                        
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPrior.R
## value   ?                                                          
## visible FALSE                                                      
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorESAndEpsFixed.R
## value   ?                                                                       
## visible FALSE                                                                   
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorNoCK.R
## value   ?                                                              
## visible FALSE                                                          
##         ./Functions/Model/ModelUtils/AltPlotSimTrainingCurves.R
## value   ?                                                      
## visible FALSE                                                  
##         ./Functions/Model/ModelUtils/aModelUtils.R
## value   ?                                         
## visible FALSE                                     
##         ./Functions/Model/ModelUtils/PlotAllWorstTestOptionsSimVsEmp.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/ModelUtils/PlotSimEmpTest.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/ModelUtils/PlotSimEmpTestNoSim.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/ModelUtils/PrepDfForModel.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/ModelUtils/PrepDfForModelPROpt.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/ModelUtils/RecodeDfs.R
## value   ?                                       
## visible FALSE                                   
##         ./Functions/Model/ModelUtils/SimplePlotSimTrainingDelay.R
## value   ?                                                        
## visible FALSE                                                    
##         ./Functions/Model/OptFxs/RunOptTrainSIT.R
## value   ?                                        
## visible FALSE                                    
##         ./Functions/Model/OptFxs/RunOptTrainSITPR.R
## value   ?                                          
## visible FALSE                                      
##         ./Functions/Model/SimFxs/RunSimTrainSIT.R
## value   ?                                        
## visible FALSE                                    
##         ./Functions/Model/SimFxs/RunSimTrainSITForPR.R
## value   ?                                             
## visible FALSE                                         
##         ./Functions/Model/TrialWiseComps/aModelFunctions.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/TrialWiseComps/CalcBayesMean.R
## value   ?                                               
## visible FALSE                                           
##         ./Functions/Model/TrialWiseComps/CalcBayesStd.R
## value   ?                                              
## visible FALSE                                          
##         ./Functions/Model/TrialWiseComps/CalcBayesVar.R
## value   ?                                              
## visible FALSE                                          
##         ./Functions/Model/TrialWiseComps/DecayChoiceKernel.R
## value   ?                                                   
## visible FALSE                                               
##         ./Functions/Model/TrialWiseComps/DecayQVals.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/TrialWiseComps/UpdateABMatrix.R
## value   ?                                                
## visible FALSE
DefPlotPars()
registerDoParallel(cores=round(detectCores()*2/3))

Load data from studies 1 and 2

s1_train <- read.csv("../data/cleaned_files/s1_train_with_delay_deident.csv")
  #read.csv("../data/cleaned_files/s1_train_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s1_sit <- read.csv("../data/cleaned_files/s1_SIT_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_train <- read.csv("../data/cleaned_files/s2_train_with_delay_deident.csv") #read.csv("../data/cleaned_files/s2_train_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_sit <- read.csv("../data/cleaned_files/s2_sit_deident_corrected_names.csv") 

Harmonize variables and create some separate vars

# Study 2 harmonize  
s2_sit[s2_sit$condition=="cogn", "condition"] <- "cognitive" 
s2_train[s2_train$trial_within_condition <= 20, "block"] <- 1
s2_train[s2_train$trial_within_condition > 20, "block"] <- 2

s1_sit$probability <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 1)))
s1_sit$valence <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 2)))

s2_sit$probability <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 1)))
s2_sit$valence <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 2)))

Define paths and functions

# MLE results  
allp <- "../model_res/opts_mle_paper_final/all/"
# Empirical Bayes results 
bp <- "../model_res/opts_mle_paper_final/best/"
# Simulations  
sp <- "../model_res/sims_clean/sims_from_mle/"
# Read model function 
rm <- function(path, model_str) read.csv(paste0(path, model_str))

Behavioral Results for Learning and Test phase

Visualization

Summarize training mean and within-subject SEM

s1_train_summs <- s1_train %>% group_by(condition, ID) %>%  summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s1_tr_summs <- Rmisc::summarySEwithin(s1_train_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s2_train_summs <- s2_train %>% group_by(condition, ID) %>%  summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s2_tr_summs <- Rmisc::summarySEwithin(s2_train_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s2_tr_summs <- Rmisc::summarySEwithin(s2_train_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s1_test_summs <- s1_sit %>% group_by(condition, ID) %>%  summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s1_te_summs <- Rmisc::summarySEwithin(s1_test_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s2_test_summs <- s2_sit %>% group_by(condition, ID) %>%  summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s2_te_summs <- Rmisc::summarySEwithin(s2_test_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s2_te_summs <- Rmisc::summarySEwithin(s2_test_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s1_tr_summs$experiment <- 1
s1_tr_summs$phase <- "Learning"

s2_tr_summs$experiment <- 2
s2_tr_summs$phase <- "Learning"

s1_te_summs$experiment <- 1
s1_te_summs$phase <- "Test"

s2_te_summs$experiment <- 2
s2_te_summs$phase <- "Test"

comb_summs <- rbind(s1_tr_summs, s2_tr_summs)
comb_summs$experiment <- factor(comb_summs$experiment)

comb_te_summs <- rbind(s1_te_summs, s2_te_summs)
comb_te_summs$experiment <- factor(comb_te_summs$experiment)
a <- ggplot(comb_summs, 
      aes(x=experiment, y=m, group=condition, fill=condition)) + 
      geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3) + 
      geom_hline(yintercept=.5, size=2) + 
      scale_fill_manual(values=c("purple", "orange")) +
      geom_errorbar(aes(x=experiment, ymin=m-se, ymax=m+se,
                        group=condition), 
                    inherit.aes=FALSE, size=1.5, width=.25, position = position_dodge(width = .3)) +
  geom_point(size=4, color="black", pch=21, position = position_dodge(width = .3)) + ga + ap + tol +
  xlab("") + ylab("Proportion correct") + ggtitle("Learning") + tp +
  ylim(c(.45, .82)) + theme(axis.text.x = element_text(size=30))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
a

b <- ggplot(comb_te_summs, 
      aes(x=experiment, y=m, group=condition, fill=condition)) + 
      geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3) + 
      geom_hline(yintercept=.5, size=2) + 
      scale_fill_manual(values=c("purple", "orange")) +
      geom_errorbar(aes(x=experiment, ymin=m-se, ymax=m+se,
                        group=condition), 
                    inherit.aes=FALSE, size=1.5, width=.25, position = position_dodge(width = .3)) +
  geom_point(size=3, color="black", pch=21, position = position_dodge(width = .3)) + ga + ap + tol + tp + 
  #lp + Used to cut legend and then turned off so size is comparable 
  xlab("") + ylab("") + ggtitle("Test") +
  ylim(c(.45, .82)) + theme(axis.text.x = element_text(size=30))
b

perf_comb <- a + b
perf_comb

#ggsave("../paper/figs/fig_parts/perf_plot.png", perf_comb, height=6, width=10, dpi=700)

Statistical analyses

Create sum contrast codes

# Specify overt as baseline factor so that negative effect corresponds to worst performance  
s1_train$cond_cc <- factor(s1_train$condition, levels=c("overt", "cognitive"))
contrasts(s1_train$cond_cc) <- c(-.5, .5)
#head(s1_train$cond_cc)

s2_train$cond_cc <- factor(s2_train$condition, levels=c("overt", "cognitive"))
contrasts(s2_train$cond_cc) <- c(-.5, .5)

s1_sit$cond_cc <- factor(s1_sit$condition, levels=c("overt", "cognitive"))
contrasts(s1_sit$cond_cc) <- c(-.5, .5)

s2_sit$cond_cc <- factor(s2_sit$condition, levels=c("overt", "cognitive"))
contrasts(s2_sit$cond_cc) <- c(-.5, .5)

Regressions of condition differences in performance

Experiment 1 — Learning and test models

summary(m1_train_full_regr <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + 
                                      (cond_cc +  scale(trial_within_condition) |ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (cond_cc +  
##     scale(trial_within_condition) | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  45054.2  45131.6 -22518.1  45036.2    39991 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.9157 -1.0686  0.4573  0.7047  1.2827 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr       
##  ID     (Intercept)                   0.4466   0.6683              
##         cond_cc1                      0.4926   0.7018   -0.19      
##         scale(trial_within_condition) 0.1156   0.3399    0.82  0.01
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.04191    0.06122  17.019  < 2e-16 ***
## cond_cc1                      -0.18054    0.06745  -2.677  0.00743 ** 
## scale(trial_within_condition)  0.36991    0.03300  11.210  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1
## cond_cc1    -0.174       
## scl(trl_w_)  0.758  0.009
summary(m1_test_full_regr <- glmer(correct ~  cond_cc + (cond_cc|ID), 
          data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
##    Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7854.4   7889.3  -3922.2   7844.4     7995 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1729  0.1319  0.2789  0.5791  1.2316 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 1.558    1.248         
##         cond_cc1    1.972    1.404    -0.17
## Number of obs: 8000, groups:  ID, 125
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.5690     0.1193  13.149  < 2e-16 ***
## cond_cc1     -0.4569     0.1510  -3.027  0.00247 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.173

Experiment 2 — Learning and test models

summary(m2_train_full_regr <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + 
                                      (cond_cc +  scale(trial_within_condition) | ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (cond_cc +  
##     scale(trial_within_condition) | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  48366.7  48445.0 -24174.4  48348.7    44151 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.0445  -1.0336   0.4405   0.6872   1.4792 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr       
##  ID     (Intercept)                   0.6677   0.8171              
##         cond_cc1                      0.5648   0.7515   -0.07      
##         scale(trial_within_condition) 0.1676   0.4094    0.87 -0.04
## Number of obs: 44160, groups:  ID, 138
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.14913    0.07101  16.184   <2e-16 ***
## cond_cc1                      -0.17464    0.06852  -2.549   0.0108 *  
## scale(trial_within_condition)  0.39251    0.03733  10.515   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1
## cond_cc1    -0.071       
## scl(trl_w_)  0.823 -0.034
summary(m2_test_full_regr <- glmer(correct ~  cond_cc + (cond_cc|ID), 
          data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
##    Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   8426.3   8461.7  -4208.1   8416.3     8827 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2283  0.1250  0.2301  0.5751  1.5541 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 1.850    1.360         
##         cond_cc1    2.467    1.571    -0.15
## Number of obs: 8832, groups:  ID, 138
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.6720     0.1235  13.542   <2e-16 ***
## cond_cc1     -0.3543     0.1586  -2.234   0.0255 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.156
car::vif(m1_train_full_regr)
##                       cond_cc scale(trial_within_condition) 
##                      1.000075                      1.000075
car::vif(m2_train_full_regr)
##                       cond_cc scale(trial_within_condition) 
##                      1.001187                      1.001187

Robustness check 1: Same models with just random intercepts instead of both random intercepts and slopes

Experiment 1

summary(m1_train_ri_only_regr <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + 
                                      (1 |ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (1 | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  46071.1  46105.5 -23031.6  46063.1    39996 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6755 -1.1011  0.5027  0.6818  1.2499 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.322    0.5675  
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    0.95856    0.05214  18.386  < 2e-16 ***
## cond_cc1                      -0.14312    0.02271  -6.302 2.93e-10 ***
## scale(trial_within_condition)  0.28711    0.01147  25.021  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1
## cond_cc1    -0.006       
## scl(trl_w_)  0.026 -0.008
summary(m1_test_ri_only_regr <- glmer(correct ~  cond_cc + (1 |ID), 
          data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (1 | ID)
##    Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   8080.6   8101.5  -4037.3   8074.6     7997 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4964  0.1516  0.3789  0.6174  1.2182 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 1.299    1.14    
## Number of obs: 8000, groups:  ID, 125
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.40802    0.10812  13.023  < 2e-16 ***
## cond_cc1    -0.31767    0.05567  -5.707 1.15e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.022

Experiment 2

summary(m2_train_ri_only_regr <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + 
                                      (1  | ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (1 | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  49635.6  49670.3 -24813.8  49627.6    44156 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4739 -1.0679  0.4853  0.6668  1.2916 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.4132   0.6428  
## Number of obs: 44160, groups:  ID, 138
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.03285    0.05596  18.458  < 2e-16 ***
## cond_cc1                      -0.15224    0.02195  -6.937 4.01e-12 ***
## scale(trial_within_condition)  0.27618    0.01108  24.922  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1
## cond_cc1    -0.007       
## scl(trl_w_)  0.024 -0.008
summary(m2_test_ri_only_regr <- glmer(correct ~  cond_cc + (1 |ID), 
          data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (1 | ID)
##    Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   8769.8   8791.0  -4381.9   8763.8     8829 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6214  0.1441  0.3650  0.6075  1.2670 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 1.577    1.256   
## Number of obs: 8832, groups:  ID, 138
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.48452    0.11314  13.121  < 2e-16 ***
## cond_cc1    -0.20840    0.05335  -3.906 9.38e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.013

Robustness check 2: Learning models with the time covariate removed

Experiments 1 and 2

summary(m1_train_no_cov_regr <- glmer(correct ~  
                                      cond_cc + 
                                      (cond_cc |ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  46182.9  46225.9 -23086.5  46172.9    39995 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3753 -1.1134  0.5036  0.6746  1.0940 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 0.3342   0.5781        
##         cond_cc1    0.4616   0.6794   -0.20
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.97023    0.05303  18.296  < 2e-16 ***
## cond_cc1    -0.17611    0.06533  -2.696  0.00702 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.189
summary(m2_train_no_cov_regr <- glmer(correct ~  
                                      cond_cc + 
                                      (cond_cc | ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  49610.0  49653.5 -24800.0  49600.0    44155 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9719 -1.0789  0.4729  0.6802  1.0720 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 0.4287   0.6548        
##         cond_cc1    0.5173   0.7192   -0.09
## Number of obs: 44160, groups:  ID, 138
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.04886    0.05698  18.407  < 2e-16 ***
## cond_cc1    -0.16968    0.06571  -2.582  0.00981 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.088
sjPlot::tab_model(m1_train_full_regr)
  correct
Predictors Odds Ratios CI p
(Intercept) 2.83 2.51 – 3.20 <0.001
cond cc1 0.83 0.73 – 0.95 0.007
trial within condition 1.45 1.36 – 1.54 <0.001
Random Effects
σ2 3.29
τ00 ID 0.45
τ11 ID.cond_cc1 0.49
τ11 ID.scale(trial_within_condition) 0.12
ρ01 -0.19
0.82
ICC 0.17
N ID 125
Observations 40000
Marginal R2 / Conditional R2 0.035 / 0.202
sjPlot::tab_model(m2_train_full_regr)
  correct
Predictors Odds Ratios CI p
(Intercept) 3.16 2.75 – 3.63 <0.001
cond cc1 0.84 0.73 – 0.96 0.011
trial within condition 1.48 1.38 – 1.59 <0.001
Random Effects
σ2 3.29
τ00 ID 0.67
τ11 ID.cond_cc1 0.56
τ11 ID.scale(trial_within_condition) 0.17
ρ01 -0.07
0.87
ICC 0.23
N ID 138
Observations 44160
Marginal R2 / Conditional R2 0.037 / 0.257
sjPlot::tab_model(m1_test_full_regr)
  correct
Predictors Odds Ratios CI p
(Intercept) 4.80 3.80 – 6.07 <0.001
cond cc1 0.63 0.47 – 0.85 0.002
Random Effects
σ2 3.29
τ00 ID 1.56
τ11 ID.cond_cc1 1.97
ρ01 ID -0.17
ICC 0.38
N ID 125
Observations 8000
Marginal R2 / Conditional R2 0.010 / 0.390
sjPlot::tab_model(m2_test_full_regr)
  correct
Predictors Odds Ratios CI p
(Intercept) 5.32 4.18 – 6.78 <0.001
cond cc1 0.70 0.51 – 0.96 0.025
Random Effects
σ2 3.29
τ00 ID 1.85
τ11 ID.cond_cc1 2.47
ρ01 ID -0.15
ICC 0.43
N ID 138
Observations 8832
Marginal R2 / Conditional R2 0.005 / 0.432

Robustness check 3: Covariate for answer type

s1_train[s1_train$response=="sum", "resp_type_cc_cog"] <- .5
s1_train[s1_train$response=="difference", "resp_type_cc_cog"] <- -.5
s1_train[s1_train$response=="top", "resp_type_cc_cog"] <- 0
s1_train[s1_train$response=="bottom", "resp_type_cc_cog"] <- 0


s1_train[s1_train$response=="top", "resp_type_cc_overt"] <- .5
s1_train[s1_train$response=="bottom", "resp_type_cc_overt"] <- -.5
s1_train[s1_train$response=="sum", "resp_type_cc_overt"] <- 0
s1_train[s1_train$response=="difference", "resp_type_cc_overt"] <- 0

s1_sit[s1_sit$resp_as_category=="sum", "resp_type_cc_cog"] <- .5
s1_sit[s1_sit$resp_as_category=="difference", "resp_type_cc_cog"] <- -.5
s1_sit[s1_sit$resp_as_category=="top", "resp_type_cc_cog"] <- 0
s1_sit[s1_sit$resp_as_category=="bottom", "resp_type_cc_cog"] <- 0

s1_sit[s1_sit$resp_as_category=="top", "resp_type_cc_overt"] <- .5
s1_sit[s1_sit$resp_as_category=="bottom", "resp_type_cc_overt"] <- -.5
s1_sit[s1_sit$resp_as_category=="sum", "resp_type_cc_overt"] <- 0
s1_sit[s1_sit$resp_as_category=="difference", "resp_type_cc_overt"] <- 0

s2_train[s2_train$response=="alphabetize", "resp_type_cc_cog"] <- .5
s2_train[s2_train$response=="rev_alphabetize", "resp_type_cc_cog"] <- -.5
s2_train[s2_train$response=="slash", "resp_type_cc_cog"] <- 0
s2_train[s2_train$response=="backslash", "resp_type_cc_cog"] <- 0

s2_train[s2_train$response=="slash", "resp_type_cc_overt"] <- .5
s2_train[s2_train$response=="backslash", "resp_type_cc_overt"] <- -.5
s2_train[s2_train$response=="alphabetize", "resp_type_cc_overt"] <- 0
s2_train[s2_train$response=="rev_alphabetize", "resp_type_cc_overt"] <- 0

s2_sit[s2_sit$resp_as_category=="alphabetize", "resp_type_cc_cog"] <- .5
s2_sit[s2_sit$resp_as_category=="rev_alphabetize", "resp_type_cc_cog"] <- -.5
s2_sit[s2_sit$resp_as_category=="slash", "resp_type_cc_cog"] <- 0
s2_sit[s2_sit$resp_as_category=="backslash", "resp_type_cc_cog"] <- 0

s2_sit[s2_sit$resp_as_category=="slash", "resp_type_cc_overt"] <- .5
s2_sit[s2_sit$resp_as_category=="backslash", "resp_type_cc_overt"] <- -.5
s2_sit[s2_sit$resp_as_category=="alphabetize", "resp_type_cc_overt"] <- 0
s2_sit[s2_sit$resp_as_category=="rev_alphabetize", "resp_type_cc_overt"] <- 0

table(s1_train$resp_type_cc_cog)
## 
##  -0.5     0   0.5 
##  9020 20000 10980
table(s1_train$resp_type_cc_overt)
## 
##  -0.5     0   0.5 
##  9451 20000 10549
table(s1_sit$resp_type_cc_cog)
## 
## -0.5    0  0.5 
## 1851 4000 2149
table(s1_sit$resp_type_cc_overt)
## 
## -0.5    0  0.5 
## 1975 4000 2025
table(s2_train$resp_type_cc_cog)
## 
##  -0.5     0   0.5 
## 10182 22080 11898
table(s2_train$resp_type_cc_overt)
## 
##  -0.5     0   0.5 
## 11186 22080 10894
table(s2_sit$resp_type_cc_cog)
## 
## -0.5    0  0.5 
## 2024 4416 2392
table(s2_sit$resp_type_cc_overt)
## 
## -0.5    0  0.5 
## 2239 4416 2177

Separately coded contrast variables for cog and overt. Random slopes of these removed for test bc singular (which makes sense bc of few observations/subj)

summary(m1_train_full_regr_at <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt +
                                      (cond_cc +  scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt |ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## correct ~ cond_cc + scale(trial_within_condition) + resp_type_cc_cog +  
##     resp_type_cc_overt + (cond_cc + scale(trial_within_condition) +  
##     resp_type_cc_cog + resp_type_cc_overt | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    45040    45212   -22500    45000    39980 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.2191 -1.0643  0.4553  0.7031  1.4380 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr                   
##  ID     (Intercept)                   0.44641  0.6681                          
##         cond_cc1                      0.49626  0.7045   -0.19                  
##         scale(trial_within_condition) 0.11598  0.3406    0.81  0.01            
##         resp_type_cc_cog              0.05454  0.2335    0.19 -0.31  0.33      
##         resp_type_cc_overt            0.07933  0.2817   -0.04  0.23  0.16  0.70
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.06182    0.06141  17.290  < 2e-16 ***
## cond_cc1                      -0.18291    0.06839  -2.674  0.00749 ** 
## scale(trial_within_condition)  0.37205    0.03308  11.248  < 2e-16 ***
## resp_type_cc_cog              -0.07227    0.04151  -1.741  0.08166 .  
## resp_type_cc_overt            -0.08409    0.04532  -1.855  0.06353 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) cnd_c1 sc(__) rsp_typ_cc_c
## cond_cc1     -0.172                           
## scl(trl_w_)   0.753  0.004                    
## rsp_typ_cc_c  0.077 -0.175  0.164             
## rsp_typ_cc_v -0.034  0.137  0.088  0.197
# Singular 
# summary(m1_test_full_regr_at <- glmer(correct ~  cond_cc + resp_type_cc_cog + resp_type_cc_overt +
#                                         (cond_cc + resp_type_cc_cog + resp_type_cc_overt |ID), 
#           data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))

summary(m1_test_full_regr_at <- glmer(correct ~  cond_cc + resp_type_cc_cog + resp_type_cc_overt +
                                        (cond_cc  |ID), 
          data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt + (cond_cc |      ID)
##    Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7857.3   7906.2  -3921.6   7843.3     7993 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0858  0.1296  0.2776  0.5823  1.2399 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 1.555    1.247         
##         cond_cc1    1.979    1.407    -0.17
## Number of obs: 8000, groups:  ID, 125
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.57061    0.11923  13.173  < 2e-16 ***
## cond_cc1           -0.45590    0.15120  -3.015  0.00257 ** 
## resp_type_cc_cog   -0.05969    0.08348  -0.715  0.47460    
## resp_type_cc_overt -0.07056    0.08936  -0.790  0.42974    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) cnd_c1 rsp_typ_cc_c
## cond_cc1     -0.175                    
## rsp_typ_cc_c -0.015 -0.021             
## rsp_typ_cc_v -0.007  0.011  0.000
summary(m2_train_full_regr_at <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt +
                                      (cond_cc +  scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt |ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## correct ~ cond_cc + scale(trial_within_condition) + resp_type_cc_cog +  
##     resp_type_cc_overt + (cond_cc + scale(trial_within_condition) +  
##     resp_type_cc_cog + resp_type_cc_overt | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  48297.6  48471.5 -24128.8  48257.6    44140 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.4618 -1.0313  0.4323  0.6859  1.4685 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr                   
##  ID     (Intercept)                   0.6757   0.8220                          
##         cond_cc1                      0.5985   0.7736   -0.07                  
##         scale(trial_within_condition) 0.1666   0.4081    0.87 -0.04            
##         resp_type_cc_cog              0.2001   0.4474    0.01 -0.26  0.01      
##         resp_type_cc_overt            0.1730   0.4160   -0.27  0.00 -0.27  0.18
## Number of obs: 44160, groups:  ID, 138
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.19401    0.07162  16.671   <2e-16 ***
## cond_cc1                      -0.16320    0.07110  -2.296   0.0217 *  
## scale(trial_within_condition)  0.39842    0.03726  10.694   <2e-16 ***
## resp_type_cc_cog              -0.12226    0.05288  -2.312   0.0208 *  
## resp_type_cc_overt            -0.00198    0.05164  -0.038   0.9694    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) cnd_c1 sc(__) rsp_typ_cc_c
## cond_cc1     -0.064                           
## scl(trl_w_)   0.821 -0.037                    
## rsp_typ_cc_c -0.001 -0.199  0.016             
## rsp_typ_cc_v -0.183  0.007 -0.170  0.087
# Singular  
# summary(m2_test_full_regr_at <- glmer(correct ~  cond_cc + resp_type_cc_cog + resp_type_cc_overt +
#                                         (cond_cc  + resp_type_cc_cog + resp_type_cc_overt  |ID), 
#           data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))

summary(m2_test_full_regr_at <- glmer(correct ~  cond_cc + resp_type_cc_cog + resp_type_cc_overt +
                                        (cond_cc  |ID), 
          data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt + (cond_cc |      ID)
##    Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   8427.9   8477.5  -4207.0   8413.9     8825 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3845  0.1211  0.2303  0.5752  1.5798 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 1.848    1.359         
##         cond_cc1    2.469    1.571    -0.16
## Number of obs: 8832, groups:  ID, 138
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.67413    0.12343  13.563   <2e-16 ***
## cond_cc1           -0.35472    0.15864  -2.236   0.0254 *  
## resp_type_cc_cog   -0.05734    0.08376  -0.685   0.4936    
## resp_type_cc_overt  0.11865    0.08580   1.383   0.1667    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) cnd_c1 rsp_typ_cc_c
## cond_cc1     -0.159                    
## rsp_typ_cc_c -0.014 -0.022             
## rsp_typ_cc_v  0.010 -0.015  0.000

Robustness check 4: paired t-test on arcsine-transformed proportion data instead of mixed-effects model

s1_summs_id <- s1_train %>% group_by(condition, ID) %>% 
   summarize(m=mean(correct)) 
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s2_summs_id <- s2_train %>% group_by(condition, ID) %>% 
   summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s1_summs_test_id <- s1_sit %>% group_by(condition, ID) %>% 
   summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s2_summs_test_id <- s2_sit %>% group_by(condition, ID) %>% 
   summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
t.test(asin(sqrt(m)) ~ condition, data=s1_summs_id, paired=TRUE)
## 
##  Paired t-test
## 
## data:  asin(sqrt(m)) by condition
## t = -2.5152, df = 124, p-value = 0.01318
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.061942537 -0.007385729
## sample estimates:
## mean of the differences 
##             -0.03466413
t.test(asin(sqrt(m)) ~ condition, data=s2_summs_id, paired=TRUE)
## 
##  Paired t-test
## 
## data:  asin(sqrt(m)) by condition
## t = -2.5139, df = 137, p-value = 0.0131
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.061281292 -0.007319414
## sample estimates:
## mean of the differences 
##             -0.03430035
t.test(asin(sqrt(m)) ~ condition, data=s1_summs_test_id, paired=TRUE)
## 
##  Paired t-test
## 
## data:  asin(sqrt(m)) by condition
## t = -2.6957, df = 124, p-value = 0.008001
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.12657517 -0.01939599
## sample estimates:
## mean of the differences 
##             -0.07298558
t.test(asin(sqrt(m)) ~ condition, data=s2_summs_test_id, paired=TRUE)
## 
##  Paired t-test
## 
## data:  asin(sqrt(m)) by condition
## t = -1.8641, df = 137, p-value = 0.06445
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.105733994  0.003119297
## sample estimates:
## mean of the differences 
##             -0.05130735

Delay effects

Create sum contrast codes for delay

# Specify overt as baseline factor so that negative effect corresponds to worst performance  
s1_train$cond_cc <- factor(s1_train$condition, levels=c("overt", "cognitive"))
contrasts(s1_train$cond_cc) <- c(-.5, .5)

s2_train$cond_cc <- factor(s2_train$condition, levels=c("overt", "cognitive"))
contrasts(s2_train$cond_cc) <- c(-.5, .5)
#head(factor(if_else(s1_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay")))
s1_train$delay_cc <- factor(if_else(s1_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay"))
contrasts(s1_train$delay_cc) <- c(-.5, .5)

s2_train$delay_cc <- factor(if_else(s2_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay"))
contrasts(s2_train$delay_cc) <- c(-.5, .5)

Building up in complexity…

Effect of delay with no moderation

summary(m1_delay <- glmer(correct ~  delay_cc + (delay_cc|ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc + (delay_cc | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  43647.9  43690.6 -21818.9  43637.9    37995 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2305 -1.1242  0.5124  0.6818  1.0042 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr
##  ID     (Intercept) 0.32784  0.5726       
##         delay_cc1   0.02447  0.1564   0.43
## Number of obs: 38000, groups:  ID, 125
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.07222    0.05325   20.14   <2e-16 ***
## delay_cc1   -0.34260    0.03240  -10.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## delay_cc1 0.043
summary(m2_delay <- glmer(correct ~  delay_cc + (delay_cc|ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc + (delay_cc | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  47084.2  47127.5 -23537.1  47074.2    41947 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2244 -1.1013  0.4908  0.6573  1.0211 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr
##  ID     (Intercept) 0.42730  0.6537       
##         delay_cc1   0.01198  0.1095   0.44
## Number of obs: 41952, groups:  ID, 138
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.13294    0.05747   19.71   <2e-16 ***
## delay_cc1   -0.30347    0.03024  -10.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## delay_cc1 -0.001

Now moderated by condition

#Singular for both 
# summary(m1_delay_cond <- glmer(correct ~  delay_cc*cond_cc + (delay_cc*cond_cc|ID), 
#           data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# 
# summary(m2_delay_cond <- glmer(correct ~  delay_cc*cond_cc + (delay_cc*cond_cc|ID), 
#           data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))

# Simplified  
summary(m1_delay_cond <- glmer(correct ~  delay_cc*cond_cc + (delay_cc + cond_cc|ID),
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc * cond_cc + (delay_cc + cond_cc | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  43090.3  43175.7 -21535.1  43070.3    37990 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0844 -1.0664  0.4738  0.6770  1.1528 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr       
##  ID     (Intercept) 0.35727  0.5977              
##         delay_cc1   0.02429  0.1559    0.37      
##         cond_cc1    0.50501  0.7106   -0.22 -0.17
## Number of obs: 38000, groups:  ID, 125
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.10896    0.05555  19.964   <2e-16 ***
## delay_cc1          -0.35536    0.03271 -10.866   <2e-16 ***
## cond_cc1           -0.17948    0.07025  -2.555   0.0106 *  
## delay_cc1:cond_cc1 -0.01904    0.05743  -0.332   0.7402    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dly_c1 cnd_c1
## delay_cc1    0.019              
## cond_cc1    -0.196 -0.055       
## dly_cc1:c_1  0.006 -0.041 -0.226
summary(m2_delay_cond <- glmer(correct ~  delay_cc*cond_cc + (delay_cc + cond_cc|ID),
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc * cond_cc + (delay_cc + cond_cc | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  46388.5  46474.9 -23184.2  46368.5    41942 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.9100 -1.0345  0.4587  0.6672  1.1512 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr       
##  ID     (Intercept) 0.45639  0.6756              
##         delay_cc1   0.01049  0.1024    0.55      
##         cond_cc1    0.56015  0.7484   -0.09 -0.21
## Number of obs: 41952, groups:  ID, 138
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.17084    0.05937  19.723   <2e-16 ***
## delay_cc1          -0.30464    0.03033 -10.043   <2e-16 ***
## cond_cc1           -0.13973    0.07013  -1.992   0.0463 *  
## delay_cc1:cond_cc1 -0.13156    0.05532  -2.378   0.0174 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dly_c1 cnd_c1
## delay_cc1    0.023              
## cond_cc1    -0.085 -0.054       
## dly_cc1:c_1  0.001 -0.025 -0.215

Model including appropriate covariates at the highest complexity that would fit (removed trial-within-cond random slope and delay*cond slope due to singular errors in studies 2 and 1 respectively)

# summary(m1_full_delay_cond <- glmer(correct ~  delay_cc*cond_cc + scale(trial_within_condition) +
#                                  (delay_cc + cond_cc + scale(trial_within_condition)|ID),
#           data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# 
# # Singular 
# summary(m2_full_delay_cond <- glmer(correct ~  delay_cc*cond_cc + scale(trial_within_condition) +
#                                  (delay_cc + cond_cc + scale(trial_within_condition)|ID),
#           data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))

# Also singular  
# summary(m1_full_delay_cond <- glmer(correct ~  delay_cc*cond_cc + scale(trial_within_condition) +
#                                  (delay_cc*cond_cc + 1|ID),
#           data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))

summary(m1_full_delay_cond <- glmer(correct ~  delay_cc*cond_cc + scale(trial_within_condition) +
                                 (delay_cc + cond_cc |ID),
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc * cond_cc + scale(trial_within_condition) +  
##     (delay_cc + cond_cc | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  42576.5  42670.5 -21277.2  42554.5    37989 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2409 -1.0187  0.4627  0.6726  1.4465 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr       
##  ID     (Intercept) 0.36544  0.6045              
##         delay_cc1   0.02941  0.1715    0.36      
##         cond_cc1    0.51948  0.7208   -0.21 -0.18
## Number of obs: 38000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.11508    0.05618  19.849   <2e-16 ***
## delay_cc1                     -0.36388    0.03352 -10.857   <2e-16 ***
## cond_cc1                      -0.18364    0.07119  -2.579   0.0099 ** 
## scale(trial_within_condition)  0.27402    0.01215  22.553   <2e-16 ***
## delay_cc1:cond_cc1            -0.01022    0.05790  -0.177   0.8599    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dly_c1 cnd_c1 sc(__)
## delay_cc1    0.029                     
## cond_cc1    -0.194 -0.064              
## scl(trl_w_)  0.018 -0.019 -0.004       
## dly_cc1:c_1  0.007 -0.041 -0.225  0.006
summary(m2_full_delay_cond <- glmer(correct ~  delay_cc*cond_cc + scale(trial_within_condition) +
                                 (delay_cc + cond_cc|ID),
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc * cond_cc + scale(trial_within_condition) +  
##     (delay_cc + cond_cc | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  45878.5  45973.6 -22928.2  45856.5    41941 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.8197 -0.9779  0.4437  0.6526  1.4299 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr       
##  ID     (Intercept) 0.46804  0.6841              
##         delay_cc1   0.01325  0.1151    0.47      
##         cond_cc1    0.57474  0.7581   -0.09 -0.18
## Number of obs: 41952, groups:  ID, 138
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.17603    0.06016  19.549   <2e-16 ***
## delay_cc1                     -0.30940    0.03083 -10.036   <2e-16 ***
## cond_cc1                      -0.14541    0.07100  -2.048   0.0405 *  
## scale(trial_within_condition)  0.26368    0.01172  22.489   <2e-16 ***
## delay_cc1:cond_cc1            -0.11653    0.05571  -2.092   0.0365 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dly_c1 cnd_c1 sc(__)
## delay_cc1    0.016                     
## cond_cc1    -0.084 -0.052              
## scl(trl_w_)  0.016 -0.013 -0.004       
## dly_cc1:c_1  0.001 -0.025 -0.214  0.010
car::vif(m1_full_delay_cond)
##                      delay_cc                       cond_cc 
##                      1.007765                      1.059549 
## scale(trial_within_condition)              delay_cc:cond_cc 
##                      1.000389                      1.056896
car::vif(m2_full_delay_cond)
##                      delay_cc                       cond_cc 
##                      1.004262                      1.051868 
## scale(trial_within_condition)              delay_cc:cond_cc 
##                      1.000268                      1.049695

Heterogeneity in condition effect

Calculate some bxal differences

Study 1

s1_perf <- s1_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_perf_s1 <- as.numeric(unlist(s1_perf[s1_perf$condition=="overt","m"]-
                    s1_perf[s1_perf$condition=="cognitive", "m"]))

# .. and test subject-level differences in performance  
s1_perf_test <- s1_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_testperf_s1 <- as.numeric(unlist(s1_perf_test[s1_perf_test$condition=="overt","m"]-
                    s1_perf_test[s1_perf_test$condition=="cognitive", "m"]))

s1_perf_overall <- s1_train %>% group_by(ID) %>% summarize(m=mean(correct))
s1_test_perf_overall <- s1_sit %>% group_by(ID) %>% summarize(m=mean(correct))

# Modest correlations across conditions 
cor.test(unlist(s1_perf[s1_perf$condition=="overt","m"]), unlist(s1_perf[s1_perf$condition=="cognitive","m"]))
## 
##  Pearson's product-moment correlation
## 
## data:  unlist(s1_perf[s1_perf$condition == "overt", "m"]) and unlist(s1_perf[s1_perf$condition == "cognitive", "m"])
## t = 5.461, df = 123, p-value = 2.507e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2885290 0.5729171
## sample estimates:
##       cor 
## 0.4417538
cor.test(unlist(s1_perf_test[s1_perf_test$condition=="overt","m"]), unlist(s1_perf_test[s1_perf_test$condition=="cognitive","m"]))
## 
##  Pearson's product-moment correlation
## 
## data:  unlist(s1_perf_test[s1_perf_test$condition == "overt", "m"]) and unlist(s1_perf_test[s1_perf_test$condition == "cognitive", "m"])
## t = 5.6112, df = 123, p-value = 1.263e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2995928 0.5809966
## sample estimates:
##       cor 
## 0.4514492

Examining if between-condition diffs are more common at the highest level of perf, because at lower ends more noise in bx. But not much relationship

cor.test(s1_perf_overall$m, o_min_c_perf_s1)
## 
##  Pearson's product-moment correlation
## 
## data:  s1_perf_overall$m and o_min_c_perf_s1
## t = 1.0074, df = 123, p-value = 0.3157
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.08652362  0.26190545
## sample estimates:
##        cor 
## 0.09045834
ComparePars(s1_perf_overall$m, o_min_c_perf_s1, use_identity_line = 0)

Study 2

s2_perf <- s2_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_perf_s2 <- as.numeric(unlist(s2_perf[s2_perf$condition=="overt","m"]-
                    s2_perf[s2_perf$condition=="cognitive", "m"]))

s2_perf_test <- s2_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_testperf_s2 <- as.numeric(unlist(s2_perf_test[s2_perf_test$condition=="overt","m"]-
                    s2_perf_test[s2_perf_test$condition=="cognitive", "m"]))


s2_perf_overall <- s2_train %>% group_by(ID) %>% summarize(m=mean(correct))
s2_test_perf_overall <- s2_sit %>% group_by(ID) %>% summarize(m=mean(correct))

cor.test(unlist(s2_perf[s2_perf$condition=="overt","m"]), unlist(s2_perf[s2_perf$condition=="cognitive","m"]))
## 
##  Pearson's product-moment correlation
## 
## data:  unlist(s2_perf[s2_perf$condition == "overt", "m"]) and unlist(s2_perf[s2_perf$condition == "cognitive", "m"])
## t = 6.0926, df = 136, p-value = 1.074e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3207662 0.5848974
## sample estimates:
##       cor 
## 0.4630508
cor.test(unlist(s2_perf_test[s2_perf_test$condition=="overt","m"]), unlist(s2_perf_test[s2_perf_test$condition=="cognitive","m"]))
## 
##  Pearson's product-moment correlation
## 
## data:  unlist(s2_perf_test[s2_perf_test$condition == "overt", "m"]) and unlist(s2_perf_test[s2_perf_test$condition == "cognitive", "m"])
## t = 4.8527, df = 136, p-value = 3.29e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2319669 0.5180282
## sample estimates:
##       cor 
## 0.3841799
cor.test(s2_perf_overall$m, o_min_c_perf_s2)
## 
##  Pearson's product-moment correlation
## 
## data:  s2_perf_overall$m and o_min_c_perf_s2
## t = -0.0051584, df = 136, p-value = 0.9959
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1675348  0.1666748
## sample estimates:
##           cor 
## -0.0004423248
ComparePars(s2_perf_overall$m, o_min_c_perf_s2, use_identity_line = 0)

Invalid relationships to perf diffs

s1_o_inv <- s1_train %>% group_by(ID) %>% summarize(o_inv=mean(overt_invalid)) 
s1_c_inv <- s1_train %>% group_by(ID) %>% summarize(c_inv=mean(cognitive_invalid)) 
assert(all(s1_o_inv$ID==s1_c_inv$ID))
s1_inv_diff <- s1_c_inv$c_inv-s1_o_inv$o_inv
cor.test(o_min_c_perf_s1, s1_inv_diff)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and s1_inv_diff
## t = 0.029401, df = 123, p-value = 0.9766
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1730371  0.1781756
## sample estimates:
##         cor 
## 0.002650977
ComparePars(o_min_c_perf_s1, s1_inv_diff, use_identity_line = 0)

s2_o_inv <- s2_train %>% group_by(ID) %>% summarize(o_inv=mean(prop_invalid_overt)) 
s2_c_inv <- s2_train %>% group_by(ID) %>% summarize(c_inv=mean(prop_invalid_cognitive)) 
assert(all(s2_o_inv$ID==s2_c_inv$ID))
s2_inv_diff <- s2_c_inv$c_inv-s2_o_inv$o_inv
cor.test(o_min_c_perf_s2, s2_inv_diff)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and s2_inv_diff
## t = 2.7091, df = 136, p-value = 0.007615
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.0614962 0.3790481
## sample estimates:
##       cor 
## 0.2262758
ComparePars(s2_inv_diff, o_min_c_perf_s2, use_identity_line = 0, xchar="cog minus overt invalid", ychar="overt minus cog correct")

Percent of pts who actually showing better performance in cognitive

Training and test study 1

# Study 1  
length(which(o_min_c_perf_s1 < 0))/length(o_min_c_perf_s1)
## [1] 0.376
length(which(o_min_c_perf_s1 == 0))/length(o_min_c_perf_s1) # 2.4% show no diff  
## [1] 0.024
length(which(o_min_c_testperf_s1 < 0))/length(o_min_c_testperf_s1)
## [1] 0.328
length(which(o_min_c_testperf_s1 == 0))/length(o_min_c_testperf_s1) 
## [1] 0.104

Training and test study 2

length(which(o_min_c_perf_s2 < 0))/length(o_min_c_perf_s2)
## [1] 0.3913043
length(which(o_min_c_perf_s2 == 0))/length(o_min_c_perf_s2) # < 1% show no diff
## [1] 0.007246377
length(which(o_min_c_testperf_s2 < 0))/length(o_min_c_testperf_s2)
## [1] 0.3550725
length(which(o_min_c_testperf_s2 == 0))/length(o_min_c_testperf_s2) 
## [1] 0.1884058

Modeling Results (Note: model numbers match those from s.R but not those used in paper…)

… because the paper numbering follows the logic of how presented there, and because we culled poor/uninformative/auxiliary models as developing the streamlined final set

Model 1 — Q learner with different decay (phi). Includes eps, CK, block explor parameter

Numbers from before bug fix on the use of prior (3/23/23) dont use: 58986 | 52624 | 57510 | 74680

m1_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior17352.csv")
m1_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior58319.csv")

m1_study1_eb <- rbind(m1_study1_eb_v1, m1_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m1_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior12123.csv")
m1_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior74336.csv")

m1_study2_eb <- rbind(m1_study2_eb_v1, m1_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

# write.csv(m1_study1_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior_merged.csv")
# write.csv(m1_study2_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior_merged.csv")
cat("\nBeta Median \n")
## 
## Beta Median
median(m1_study1_eb$beta)
## [1] 2.330957
cat("SD \n")
## SD
sd(m1_study1_eb$beta)
## [1] 5.393476
cat("\nRL LR \n")
## 
## RL LR
median(m1_study1_eb$q_LR)
## [1] 0.5916992
cat("SD \n")
## SD
sd(m1_study1_eb$q_LR)
## [1] 0.2752539
cat("\nphi cog \n")
## 
## phi cog
median(m1_study1_eb$cog_phi)
## [1] 0.1593268
cat("SD \n")
## SD
sd(m1_study1_eb$cog_phi)
## [1] 0.2250586
cat("\nphi overt \n")
## 
## phi overt
median(m1_study1_eb$overt_phi)
## [1] 0.1328623
cat("SD \n")
## SD
sd(m1_study1_eb$overt_phi)
## [1] 0.169769
cat("\nepsilon \n")
## 
## epsilon
median(m1_study1_eb$epsilon)
## [1] 0.005554817
cat("SD \n")
## SD
sd(m1_study1_eb$epsilon)
## [1] 0.02699492
cat("\nChoice LR \n")
## 
## Choice LR
median(m1_study1_eb$choice_LR)
## [1] 0.1198229
cat("SD \n")
## SD
sd(m1_study1_eb$choice_LR)
## [1] 0.1222895
cat("\nES \n")
## 
## ES
median(m1_study1_eb$explor_scalar)
## [1] 0.3783385
cat("SD \n")
## SD
sd(m1_study1_eb$explor_scalar)
## [1] 0.2089571
cat("\nBeta Median \n")
## 
## Beta Median
median(m1_study2_eb$beta)
## [1] 2.692798
cat("SD \n")
## SD
sd(m1_study2_eb$beta)
## [1] 4.910113
cat("\nRL LR \n")
## 
## RL LR
median(m1_study2_eb$q_LR)
## [1] 0.5543408
cat("SD \n")
## SD
sd(m1_study2_eb$q_LR)
## [1] 0.2806
cat("\nphi cog \n")
## 
## phi cog
median(m1_study2_eb$cog_phi)
## [1] 0.2033696
cat("SD \n")
## SD
sd(m1_study2_eb$cog_phi)
## [1] 0.2858507
cat("\nphi overt \n")
## 
## phi overt
median(m1_study2_eb$overt_phi)
## [1] 0.157541
cat("SD \n")
## SD
sd(m1_study2_eb$overt_phi)
## [1] 0.2223507
cat("\nepsilon \n")
## 
## epsilon
median(m1_study2_eb$epsilon)
## [1] 7.844961e-07
cat("SD \n")
## SD
sd(m1_study2_eb$epsilon)
## [1] 0.02668559
cat("\nChoice LR \n")
## 
## Choice LR
median(m1_study2_eb$choice_LR)
## [1] 0.09913484
cat("SD \n")
## SD
sd(m1_study2_eb$choice_LR)
## [1] 0.128331
cat("\nES \n")
## 
## ES
median(m1_study2_eb$explor_scalar)
## [1] 0.2382697
cat("SD \n")
## SD
sd(m1_study2_eb$explor_scalar)
## [1] 0.1951073
ComparePars(m1_study1_eb$cog_phi, m1_study1_eb$overt_phi, 
            "Phi", "Cog", "Overt")

ComparePars(m1_study2_eb$cog_phi, m1_study2_eb$overt_phi, 
            "Phi", "Cog", "Overt")

Relationships at both train and test in studies 1 and 2

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m1_study1_eb$ID)

o_min_c_phi_s1_m1 <- m1_study1_eb$overt_phi - m1_study1_eb$cog_phi

cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m1)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_phi_s1_m1
## t = -6.9629, df = 123, p-value = 1.765e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6469195 -0.3927853
## sample estimates:
##        cor 
## -0.5317171
cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m1)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_phi_s1_m1
## t = -6.0126, df = 123, p-value = 1.926e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6018399 -0.3284889
## sample estimates:
##        cor 
## -0.4766035
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m1_study2_eb$ID)

o_min_c_phi_s2_m1 <- m1_study2_eb$overt_phi - m1_study2_eb$cog_phi

cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m1)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_phi_s2_m1
## t = -4.7721, df = 136, p-value = 4.642e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5133414 -0.2259169
## sample estimates:
##        cor 
## -0.3787242
cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m1)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_phi_s2_m1
## t = -7.1401, df = 136, p-value = 5.093e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6339502 -0.3889985
## sample estimates:
##        cor 
## -0.5221611

Descriptives and chi-square tests

Higher median decay in cog in both

# hist(sqrt(m1_study1_eb$cog_phi), breaks=100)
# hist(m1_study1_eb$cog_phi, breaks=100)
# hist(m1_study1_eb$overt_phi, breaks=100)
cat("\n Study 1 cog\n")
## 
##  Study 1 cog
median(m1_study1_eb$cog_phi)
## [1] 0.1593268
cat("\n Study 1 overt\n")
## 
##  Study 1 overt
median(m1_study1_eb$overt_phi)
## [1] 0.1328623
m1_s1_phi <- rbind(
  data.frame("phi"=m1_study1_eb$overt_phi, "cond"="overt"),
  data.frame("phi"=m1_study1_eb$cog_phi, "cond"="cognitive")
)

m1_s2_phi <- rbind(
  data.frame("phi"=m1_study2_eb$overt_phi, "cond"="overt"),
  data.frame("phi"=m1_study2_eb$cog_phi, "cond"="cognitive")
)
m1_s1_phi$Experiment <- "Experiment 1"
m1_s2_phi$Experiment <- "Experiment 2"

m1_phi_ID_comb <- rbind(m1_s1_phi, m1_s2_phi)

phi_medians <- rbind(
  data.frame(m1_s1_phi %>% group_by(cond) %>% summarize(m=median(phi)), "Experiment"="Experiment 1"),
  data.frame(m1_s2_phi %>% group_by(cond) %>% summarize(m=median(phi)), "Experiment"="Experiment 2")
)

Difference higher for means but not as meaningful bc of skew, so plotting median

m1_s1_phi %>% group_by(cond) %>% summarize(m=mean(phi))
## # A tibble: 2 × 2
##   cond          m
##   <chr>     <dbl>
## 1 cognitive 0.235
## 2 overt     0.176
m1_s2_phi %>% group_by(cond) %>% summarize(m=mean(phi))
## # A tibble: 2 × 2
##   cond          m
##   <chr>     <dbl>
## 1 cognitive 0.310
## 2 overt     0.209
phi_plot <-ggplot(phi_medians, aes(x=cond, y=m, color=cond, group=cond)) +
  
  geom_jitter(data = m1_phi_ID_comb, 
              aes(x=cond, y=phi, group=cond), 
              color="black", fill="gray57", size=2, pch=21, width=.15) + 
  geom_bar(stat="identity", fill="white", size=3, alpha=.2) +
  facet_wrap(~ Experiment)  + ga + ap + ft + lp + ylab(TeX('$\\phi') ) + xlab("") + 
  scale_color_manual(values=c("purple", "orange")) + ylim(-.03, 1.1) + tol
phi_plot

#ggsave("../paper/figs/fig_parts/phi_plot.png", phi_plot, height=5, width=10, dpi=700)

Percent higher phi in studies 1 and 2

c_min_o_phi_s1_m1 <- m1_study1_eb$cog_phi - m1_study1_eb$overt_phi
c_min_o_phi_s2_m1 <- m1_study2_eb$cog_phi - m1_study2_eb$overt_phi
table((c_min_o_phi_s1_m1 > 0)*1)[2]/sum(table((c_min_o_phi_s1_m1 > 0)*1))
##     1 
## 0.576
table((c_min_o_phi_s2_m1 > 0)*1)[2]/sum(table((c_min_o_phi_s2_m1 > 0)*1))
##         1 
## 0.6304348
chisq.test(table((c_min_o_phi_s1_m1 > 0)*1))
## 
##  Chi-squared test for given probabilities
## 
## data:  table((c_min_o_phi_s1_m1 > 0) * 1)
## X-squared = 2.888, df = 1, p-value = 0.08924
chisq.test(table((c_min_o_phi_s2_m1 > 0)*1))
## 
##  Chi-squared test for given probabilities
## 
## data:  table((c_min_o_phi_s2_m1 > 0) * 1)
## X-squared = 9.3913, df = 1, p-value = 0.00218
m1_s1_plot_df <- data.frame("train_oc_diff"=o_min_c_perf_s1, "test_oc_diff"=o_min_c_testperf_s1, "phi_co_diff"=c_min_o_phi_s1_m1)
m1_s2_plot_df <- data.frame("train_oc_diff"=o_min_c_perf_s2, "test_oc_diff"=o_min_c_testperf_s2, "phi_co_diff"=c_min_o_phi_s2_m1)

Training plots

a <- ggplot(m1_s1_plot_df, aes(x=phi_co_diff, y=train_oc_diff)) + 
  geom_smooth(method="lm", se=FALSE, size=3, color="black") +
  geom_point(pch=21, fill="gray57", size=6) + 
  stat_cor(method="pearson", size=6, label.y=.5) + 
  ylab("Proportion correct: \nOvert minus cognitive") + 
  xlab(TeX("")) +
  ga + ap + 
  xlim(-.9, 1.05) + ylim(-.7, .65)
a
## `geom_smooth()` using formula = 'y ~ x'

b <- ggplot(m1_s2_plot_df, aes(x=phi_co_diff, y=train_oc_diff)) + 
  geom_smooth(method="lm", se=FALSE, size=3, color="black") +
  geom_point(pch=21, fill="gray57", size=6) + 
  stat_cor(method="pearson", size=6, label.y=.5) + 
  ylab("") + 
  xlab(TeX("")) +
  ga + ap + 
  xlim(-.9, 1.05) + ylim(-.7, .65)
b
## `geom_smooth()` using formula = 'y ~ x'

Test plots

c <- ggplot(m1_s1_plot_df, aes(x=phi_co_diff, y=test_oc_diff)) + 
  geom_smooth(method="lm", se=FALSE, size=3, color="gray57") +
  geom_point(pch=21, fill="gray79", size=6) + 
  stat_cor(method="pearson", size=6, label.y=.5) + 
  ylab("Proportion correct: \nOvert minus cognitive") + 
  xlab(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
  ga + ap + 
  xlim(-.9, 1.05) + ylim(-.7, .65)

d <- ggplot(m1_s2_plot_df, aes(x=phi_co_diff, y=test_oc_diff)) + 
  geom_smooth(method="lm", se=FALSE, size=3, color="gray57") +
  geom_point(pch=21, fill="gray79", size=6) + 
  stat_cor(method="pearson", size=6, label.y=.5) + 
  ylab("") + 
  xlab(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
  ga + ap + 
  xlim(-.9, 1.05) + ylim(-.7, .65)
all_vs <- (a + b) / (c + d)
all_vs
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

#ggsave("../paper/figs/fig_parts/phi_vs_propdiff.png", all_vs, height=8, width=14, dpi=700)

Correlation with decay effect

cor.test(m1_study1_eb_v1$cog_phi+m1_study1_eb_v1$overt_phi, ranef(m1_full_delay_cond)$ID$delay_cc1)
## 
##  Pearson's product-moment correlation
## 
## data:  m1_study1_eb_v1$cog_phi + m1_study1_eb_v1$overt_phi and ranef(m1_full_delay_cond)$ID$delay_cc1
## t = -3.7766, df = 123, p-value = 0.0002463
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4712759 -0.1555427
## sample estimates:
##        cor 
## -0.3223455
cor.test(m1_study2_eb_v1$cog_phi+m1_study2_eb_v1$overt_phi, ranef(m2_full_delay_cond)$ID$delay_cc1)
## 
##  Pearson's product-moment correlation
## 
## data:  m1_study2_eb_v1$cog_phi + m1_study2_eb_v1$overt_phi and ranef(m2_full_delay_cond)$ID$delay_cc1
## t = -5.6966, df = 136, p-value = 7.244e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5646095 -0.2933251
## sample estimates:
##       cor 
## -0.438916
ComparePars((m1_study1_eb_v1$cog_phi+m1_study1_eb_v1$overt_phi), ranef(m1_full_delay_cond)$ID$delay_cc1, use_identity_line = 0)

ComparePars((m1_study2_eb_v1$cog_phi+m1_study2_eb_v1$overt_phi), ranef(m2_full_delay_cond)$ID$delay_cc1, use_identity_line = 0)

Model validation, main text parts

6/5/23 — reran sims to make sure completely updated/bug fixed version of par results used

s1_train_sim_m1_eb <- 
  rm(sp, 
     "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDiffDecayToPessPrior57088.csv")

s2_train_sim_m1_eb <- 
  rm(sp, 
    "SIM_EMPIRICAL_BAYES_study_2_train_SIT__train_RunMQLearnerDiffDecayToPessPrior28414.csv")

s1_test_sim_m1_eb <- 
  rm(sp, 
     "SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior57088.csv")

s2_test_sim_m1_eb <- 
  rm(sp, 
    "SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior28414.csv")
a <- AltOneSimEmpPlotRep(AltPlotSimTrainingCurves(emp_df=s1_train, s1_train_sim_m1_eb))
b <- AltOneSimEmpPlotRep(AltPlotSimTrainingCurves(emp_df=s2_train, s2_train_sim_m1_eb))
c <- PlotSimEmpTest(emp_df=s1_sit, sim_df=s1_test_sim_m1_eb)
d <- PlotSimEmpTest(emp_df=s2_sit, sim_df=s2_test_sim_m1_eb)
a

b

sim_test_comb <- c + d 
#(temporarily flipped tol -> lp to get the legend for paper)  
sim_test_comb

#ggsave("../paper/figs/fig_parts/train_sim_experiment1_plot.png", a, height=7, width=12, dpi=700)
ggsave("../paper/figs/fig_parts/test_sim_experiment1_plot.png", c, height=5, width=11, dpi=700)
ggsave("../paper/figs/fig_parts/test_sim_experiment2_plot.png", d, height=5, width=11, dpi=700)

For supplemental because of space constraints

#ggsave("../paper/figs/fig_parts/train_sim_experiment2_plot.png", b, height=7, width=12, dpi=700)

Model 3 — Same as m1 but not letting decay vary

Prior to vactoin — bug fixed files swapped in (see to_do notes)

m3_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior35433.csv")
m3_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior22212.csv")

m3_study1_eb <- rbind(m3_study1_eb_v1, m3_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m3_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior29818.csv")
m3_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior56742.csv")

m3_study2_eb <- rbind(m3_study2_eb_v1,m3_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

# write.csv(m3_study1_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior_merged.csv")
# write.csv(m3_study2_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior_merged.csv")

Worse in AIC

ComparePars(m3_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m3", "m1")

ComparePars(m3_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m3", "m1")

sum(m3_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 134.026
sum(m3_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 147.2234

Sanity check that the amount of model improvement when letting decay varies correlates with the diff in phi between conditions

assert(m1_study1_eb$ID==m3_study1_eb$ID)
assert(m1_study2_eb$ID==m3_study2_eb$ID)
cor.test(m3_study1_eb$AIC-m1_study1_eb$AIC, abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi), method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  m3_study1_eb$AIC - m1_study1_eb$AIC and abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi)
## S = 107038, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6711582
cor.test(m3_study2_eb$AIC-m1_study2_eb$AIC, abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi), method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  m3_study2_eb$AIC - m1_study2_eb$AIC and abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi)
## S = 163526, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6266436

(Function prints Pearson’s r but more appropriate Spearman’s \(\rho\) reported one cell above)

ComparePars(m1_study1_eb$AIC-m3_study1_eb$AIC, abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi), "", "AIC change", "Overt min Cog phi", use_identity_line = 0)

ComparePars(m1_study2_eb$AIC-m3_study2_eb$AIC, abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi), "", "AIC change", "Overt min Cog phi", use_identity_line = 0)

Model 4 — Same as m3 but to 0 inits

m4_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits47402.csv")
m4_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits84391.csv")

m4_study1_eb <- rbind(m4_study1_eb_v1, m4_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m4_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits41957.csv")
m4_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits37660.csv")

m4_study2_eb <- rbind(m4_study2_eb_v1, 
                      m4_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

# write.csv(m4_study1_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits_merged.csv")
# write.csv(m4_study2_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits_merged.csv")

Worse in AIC than pessimistic priors

ComparePars(m4_study1_eb$AIC, m3_study1_eb$AIC, "AIC", "m3", "m1")

ComparePars(m4_study2_eb$AIC, m3_study2_eb$AIC, "AIC", "m3", "m1")

sum(m4_study1_eb$AIC-m3_study1_eb$AIC)
## [1] 120.8681
sum(m4_study2_eb$AIC-m3_study2_eb$AIC)
## [1] 446.3359

Model 5 — Same as m1 but lets beta vary instead of decay

4/6 — bug fixed

m5_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta67330.csv")
m5_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta37484.csv")

m5_study1_eb <- rbind(m5_study1_eb_v1, m5_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m5_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta76602.csv")
m5_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta11129.csv")

m5_study2_eb <- rbind(m5_study2_eb_v1, m5_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))

Worse in AIC than diff decay

ComparePars(m5_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m5", "m1")

ComparePars(m5_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m5", "m1")

sum(m5_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 672.2797
sum(m5_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 857.6129

Model 6 — Same as m1 but lets learning rate vary instead of decay

4/6 — bug fix files swapped in

m6_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR12095.csv")
m6_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR13604.csv")

m6_study1_eb <- rbind(m6_study1_eb_v1, m6_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m6_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR69123.csv")
m6_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR43205.csv")

m6_study2_eb <- rbind(m6_study2_eb_v1, m6_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))

Worse in AIC than decay allowed to vary by condition

ComparePars(m6_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m6", "m1")

ComparePars(m6_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m6", "m1")

sum(m6_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 131.9173
sum(m6_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 173.9311

Higher overt learning rate in both

median(m6_study1_eb$cog_q_LR)
## [1] 0.5048065
median(m6_study1_eb$overt_q_LR)
## [1] 0.556566
median(m6_study2_eb$cog_q_LR)
## [1] 0.4577678
median(m6_study1_eb$overt_q_LR)
## [1] 0.556566

Learning rate correlated across conditions

ComparePars(m6_study1_eb$cog_q_LR, m6_study1_eb$overt_q_LR, 
            "q_LR", "Cog", "Overt")

ComparePars(m6_study2_eb$cog_q_LR, m6_study2_eb$overt_q_LR, 
            "q_LR", "Cog", "Overt")

Correlated with perf diffs but not as strongly as phi

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m6_study1_eb$ID)

o_min_c_q_LR_s1_m6 <- m6_study1_eb$overt_q_LR - m6_study1_eb$cog_q_LR

cor.test(o_min_c_perf_s1, o_min_c_q_LR_s1_m6)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_q_LR_s1_m6
## t = 4.361, df = 123, p-value = 2.708e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2034049 0.5088491
## sample estimates:
##       cor 
## 0.3659413
cor.test(o_min_c_testperf_s1, o_min_c_q_LR_s1_m6)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_q_LR_s1_m6
## t = 2.6939, df = 123, p-value = 0.008049
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.06303996 0.39525878
## sample estimates:
##       cor 
## 0.2360345
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m6_study2_eb$ID)

o_min_c_q_LR_s2_m6 <- m6_study2_eb$overt_q_LR - m6_study2_eb$cog_q_LR

cor.test(o_min_c_perf_s2, o_min_c_q_LR_s2_m6)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_q_LR_s2_m6
## t = 3.83, df = 136, p-value = 0.0001949
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1528938 0.4553870
## sample estimates:
##       cor 
## 0.3120266
cor.test(o_min_c_testperf_s2, o_min_c_q_LR_s2_m6)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_q_LR_s2_m6
## t = 5.4029, df = 136, p-value = 2.849e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2724024 0.5489175
## sample estimates:
##      cor 
## 0.420372
chisq.test(table((o_min_c_q_LR_s1_m6 > 0)*1))
## 
##  Chi-squared test for given probabilities
## 
## data:  table((o_min_c_q_LR_s1_m6 > 0) * 1)
## X-squared = 6.728, df = 1, p-value = 0.009491
chisq.test(table((o_min_c_q_LR_s2_m6 > 0)*1))
## 
##  Chi-squared test for given probabilities
## 
## data:  table((o_min_c_q_LR_s2_m6 > 0) * 1)
## X-squared = 1.8551, df = 1, p-value = 0.1732

Model 11 — Simple q-learner model

m11_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearner24681.csv")
m11_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearner77403.csv")

m11_study1_eb <- rbind(m11_study1_eb_v1,m11_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))


m11_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearner18790.csv")
m11_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearner26119.csv")

m11_study2_eb <- rbind(m11_study2_eb_v1, m11_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))

Much worse in AIC to not include decay

ComparePars(m11_study1_eb$AIC, m3_study1_eb$AIC, "AIC", "m11", "m3")

ComparePars(m11_study2_eb$AIC, m3_study2_eb$AIC, "AIC", "m11", "m3")

sum(m11_study1_eb$AIC-m3_study1_eb$AIC)
## [1] 918.2861
sum(m11_study2_eb$AIC-m3_study2_eb$AIC)
## [1] 1225.945

Model 12 — Simple Bayes model

m12_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner13415.csv")
m12_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner42048.csv")

m12_study1_eb <- rbind(m12_study1_eb_v1,m12_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m12_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner22977.csv")
m12_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner23536.csv")

m12_study2_eb <- rbind(m12_study2_eb_v1, m12_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))

In turn much worse than the simple Q-learning model

ComparePars(m12_study1_eb$AIC, m11_study1_eb$AIC, "AIC", "m12", "m11")

ComparePars(m12_study2_eb$AIC, m11_study2_eb$AIC, "AIC", "m12", "m11")

sum(m12_study1_eb$AIC-m11_study1_eb$AIC)
## [1] 622.5435
sum(m12_study2_eb$AIC-m11_study2_eb$AIC)
## [1] 1012.796

Model 13 — M12 but allowing beta to vary

m13_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta60557.csv")
m13_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta85855.csv")

m13_study1_eb <- rbind(m13_study1_eb_v1, m13_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m13_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta29859.csv")
m13_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta29859.csv")

m13_study2_eb <- rbind(m13_study2_eb_v1, m13_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

Way worse than m1

ComparePars(m13_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m13", "m1")

ComparePars(m13_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m13", "m1")

sum(m13_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 2081.633
sum(m13_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 2997.104

Higher median beta in overt

median(m13_study1_eb$cog_beta)
## [1] 1.790819
median(m13_study1_eb$overt_beta)
## [1] 1.876175
median(m13_study2_eb$cog_beta)
## [1] 1.893116
median(m13_study2_eb$overt_beta)
## [1] 2.099269

Modest correlations between beta and perf

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m13_study1_eb$ID)

o_min_c_beta_s1_m13 <- m13_study1_eb$overt_beta - m13_study1_eb$cog_beta

cor.test(o_min_c_perf_s1, o_min_c_beta_s1_m13)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_beta_s1_m13
## t = 3.5966, df = 123, p-value = 0.0004653
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1404791 0.4592086
## sample estimates:
##       cor 
## 0.3084768
cor.test(o_min_c_testperf_s1, o_min_c_beta_s1_m13)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_beta_s1_m13
## t = 2.7331, df = 123, p-value = 0.007199
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.06646281 0.39815509
## sample estimates:
##       cor 
## 0.2392776
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m13_study2_eb$ID)

o_min_c_beta_s2_m13 <- m13_study2_eb$overt_beta - m13_study2_eb$cog_beta

cor.test(o_min_c_perf_s2, o_min_c_beta_s2_m13)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_beta_s2_m13
## t = 2.8955, df = 136, p-value = 0.004412
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.07696563 0.39227998
## sample estimates:
##       cor 
## 0.2409713
cor.test(o_min_c_testperf_s2, o_min_c_beta_s2_m13)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_beta_s2_m13
## t = 1.9265, df = 136, p-value = 0.05613
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.004233115  0.321339865
## sample estimates:
##      cor 
## 0.162987

Model 14 — Hybrid model — integrative plus bayes — bayes contribution varies by cond

Alt refers to this being the simpler version with uncertainty factored in at the action (arm) level rather than a weighted uncertainty computed — this and the more complex version were v similar in terms of fit

m14_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt39650.csv")
m14_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt71571.csv")

m14_study1_eb <- rbind(m14_study1_eb_v1, m14_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m14_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt38226.csv")
m14_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt28286.csv")

m14_study2_eb <- rbind(m14_study2_eb_v1, m14_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

Original version with the more complicated weighting scheme

m14_study1_eb_v1_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes25885.csv")
m14_study1_eb_v2_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes65430.csv")

m14_study1_eb_old <- rbind(m14_study1_eb_v1_old, m14_study1_eb_v2_old) %>% 
  group_by(ID) %>% slice(which.min(nll))

m14_study2_eb_v1_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes30853.csv")
m14_study2_eb_v2_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes58630.csv")

m14_study2_eb_old <- rbind(m14_study2_eb_v1_old, m14_study2_eb_v2_old) %>% 
  group_by(ID) %>% slice(which.min(nll))

As expected, these are very similar, with the new version a tad better

ComparePars(m14_study1_eb$AIC, m14_study1_eb_old$AIC, "AIC", "m14", "m14 old")

sum(m14_study1_eb$AIC-m14_study1_eb_old$AIC)
## [1] -18.67696
ComparePars(m14_study2_eb$AIC, m14_study2_eb_old$AIC, "AIC", "m14", "m14 old")

sum(m14_study2_eb$AIC-m14_study2_eb_old$AIC)
## [1] -43.19563

Worse AIC (and nll) compared to m1 (not a nested model bc doesn’t have diff decay)

ComparePars(m14_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m14", "m1")

ComparePars(m14_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m14", "m1")

sum(m14_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 448.2412
sum(m14_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 486.6508
sum(m14_study1_eb$nll-m1_study1_eb$nll)
## [1] 99.12058
sum(m14_study2_eb$nll-m1_study2_eb$nll)
## [1] 105.3254

Higher median rho in overt

median(m14_study1_eb$rho_cog)
## [1] 0.1651074
median(m14_study1_eb$rho_overt)
## [1] 0.2611568
median(m14_study2_eb$rho_cog)
## [1] 0.1066345
median(m14_study2_eb$rho_overt)
## [1] 0.2367347

Correlations w perf but not as strongly as phi

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m14_study1_eb$ID)

o_min_c_beta_s1_m14 <- m14_study1_eb$rho_overt - m14_study1_eb$rho_cog

cor.test(o_min_c_perf_s1, o_min_c_beta_s1_m14)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_beta_s1_m14
## t = 4.2195, df = 123, p-value = 4.71e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1919767 0.4999821
## sample estimates:
##       cor 
## 0.3555962
cor.test(o_min_c_testperf_s1, o_min_c_beta_s1_m14)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_beta_s1_m14
## t = 3.3007, df = 123, p-value = 0.001262
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1154251 0.4388741
## sample estimates:
##       cor 
## 0.2852507
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m14_study2_eb$ID)

o_min_c_beta_s2_m14 <- m14_study2_eb$rho_overt - m14_study2_eb$rho_cog

cor.test(o_min_c_perf_s2, o_min_c_beta_s2_m14)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_beta_s2_m14
## t = 5.1757, df = 136, p-value = 7.977e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2558970 0.5363993
## sample estimates:
##       cor 
## 0.4056554
cor.test(o_min_c_testperf_s2, o_min_c_beta_s2_m14)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_beta_s2_m14
## t = 3.4613, df = 136, p-value = 0.0007187
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1232954 0.4311429
## sample estimates:
##       cor 
## 0.2845379

Model 15 — Hybrid model — integrative plus bayes — decay contribution varies by cond

Alternate version with the simpler action-wise weighting (see m14 description)

m15_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt33317.csv")
m15_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt35801.csv")

m15_study1_eb <- rbind(m15_study1_eb_v1, m15_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))


m15_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt64665.csv")
m15_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt54458.csv")

m15_study2_eb <- rbind(m15_study2_eb_v1, m15_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

Original version with the more complicated weighting scheme

bug numbers: 75340 81994 63931 40889

m15_study1_eb_v1_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes50404.csv")
m15_study1_eb_v2_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes56900.csv")

m15_study1_eb_old <- rbind(m15_study1_eb_v1_old, m15_study1_eb_v2_old) %>% 
  group_by(ID) %>% slice(which.min(nll))


m15_study2_eb_v1_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes68726.csv")
m15_study2_eb_v2_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes30103.csv")

m15_study2_eb_old <- rbind(m15_study2_eb_v1_old, m15_study2_eb_v2_old) %>% 
  group_by(ID) %>% slice(which.min(nll))

As expected, again very similar — so using new/simpler one going forward

ComparePars(m15_study1_eb$AIC, m15_study1_eb_old$AIC, "AIC", "m15", "m15 old")

sum(m15_study1_eb$AIC-m15_study1_eb_old$AIC)
## [1] 11.82675
ComparePars(m15_study2_eb$AIC, m15_study2_eb_old$AIC, "AIC", "m15", "m15 old")

sum(m15_study2_eb$AIC-m15_study2_eb_old$AIC)
## [1] -53.60774

Worse AIC compared to m1

ComparePars(m15_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m15", "m1")

ComparePars(m15_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m15", "m1")

sum(m15_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 169.3991
sum(m15_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 127.4026
# sum(m15_study1_eb$nll-m1_study1_eb$nll)
# sum(m15_study2_eb$nll-m1_study2_eb$nll)

Higher median decay in cog in both

median(m15_study1_eb$cog_phi)
## [1] 0.2204569
median(m15_study1_eb$overt_phi)
## [1] 0.1742236
median(m15_study2_eb$cog_phi)
## [1] 0.3308409
median(m15_study1_eb$overt_phi)
## [1] 0.1742236
ComparePars(m15_study1_eb$cog_phi, m15_study1_eb$overt_phi, 
            "Phi", "Cog", "Overt")

ComparePars(m15_study2_eb$cog_phi, m15_study2_eb$overt_phi, 
            "Phi", "Cog", "Overt")

Similar perf correlations as m1

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m15_study1_eb$ID)

o_min_c_phi_s1_m15 <- m15_study1_eb$overt_phi - m15_study1_eb$cog_phi

cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m15)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_phi_s1_m15
## t = -7.1705, df = 123, p-value = 6.097e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6560039 -0.4060503
## sample estimates:
##       cor 
## -0.542943
cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m15)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_phi_s1_m15
## t = -6.1734, df = 123, p-value = 8.902e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6098839 -0.3397791
## sample estimates:
##        cor 
## -0.4863663
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m15_study2_eb$ID)

o_min_c_phi_s2_m15 <- m15_study2_eb$overt_phi - m15_study2_eb$cog_phi

cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m15)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_phi_s2_m15
## t = -3.8478, df = 136, p-value = 0.0001825
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4565342 -0.1543079
## sample estimates:
##        cor 
## -0.3133333
cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m15)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_phi_s2_m15
## t = -6.5236, df = 136, p-value = 1.247e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6058791 -0.3496199
## sample estimates:
##        cor 
## -0.4882024

Model 27 — Fixing both ES and epsilon at the group level

m27_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed46645.csv")
m27_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed33416.csv")

m27_study1_eb <- rbind(m27_study1_eb_v1, m27_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m27_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed20328.csv")
m27_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed44601.csv")

m27_study2_eb <- rbind(m27_study2_eb_v1, m27_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

# write.csv(m27_study1_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed_merged.csv")
# write.csv(m27_study2_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed_merged.csv")

Does worsen NLL substantially

ComparePars(m27_study1_eb$nll, m1_study1_eb$nll, "nll", "m27", "m1")

ComparePars(m27_study2_eb$nll, m1_study2_eb$nll, "nll", "m27", "m1")

sum(m27_study1_eb$nll-m1_study1_eb$nll)
## [1] 259.0197
sum(m27_study2_eb$nll-m1_study2_eb$nll)
## [1] 367.588

Still higher median decay in cog in both

median(m27_study1_eb$cog_phi)
## [1] 0.1637271
median(m27_study1_eb$overt_phi)
## [1] 0.1337206
median(m27_study2_eb$cog_phi)
## [1] 0.2304866
median(m27_study1_eb$overt_phi)
## [1] 0.1337206

Strong relationships at both train and test in studies 1 and 2

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m27_study1_eb$ID)

o_min_c_phi_s1_m27 <- m27_study1_eb$overt_phi - m27_study1_eb$cog_phi

cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m27)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_phi_s1_m27
## t = -6.9726, df = 123, p-value = 1.68e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6473518 -0.3934142
## sample estimates:
##        cor 
## -0.5322504
cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m27)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_phi_s1_m27
## t = -6.4826, df = 123, p-value = 1.966e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6248649 -0.3610156
## sample estimates:
##       cor 
## -0.504631
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m27_study2_eb$ID)

o_min_c_phi_s2_m27 <- m27_study2_eb$overt_phi - m27_study2_eb$cog_phi

cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m27)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_phi_s2_m27
## t = -4.8922, df = 136, p-value = 2.775e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5203133 -0.2349247
## sample estimates:
##        cor 
## -0.3868433
cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m27)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_phi_s2_m27
## t = -7.114, df = 136, p-value = 5.847e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6328069 -0.3873771
## sample estimates:
##       cor 
## -0.520771

59 and 62% of pts show a higher phi in overt

table((o_min_c_phi_s1_m27 > 0)*1)
## 
##  0  1 
## 73 52
table((o_min_c_phi_s2_m27 > 0)*1)
## 
##  0  1 
## 86 52
table((o_min_c_phi_s1_m27 > 0))[1]/sum(table((o_min_c_phi_s1_m27 > 0)))
## FALSE 
## 0.584
table((o_min_c_phi_s2_m27 > 0))[1]/sum(table((o_min_c_phi_s2_m27 > 0)))
##     FALSE 
## 0.6231884
chisq.test(table(if_else((o_min_c_phi_s1_m27 > 0) == TRUE, 1, 0)))
## 
##  Chi-squared test for given probabilities
## 
## data:  table(if_else((o_min_c_phi_s1_m27 > 0) == TRUE, 1, 0))
## X-squared = 3.528, df = 1, p-value = 0.06034
chisq.test(table(if_else((o_min_c_phi_s2_m27 > 0) == TRUE, 1, 0)))
## 
##  Chi-squared test for given probabilities
## 
## data:  table(if_else((o_min_c_phi_s2_m27 > 0) == TRUE, 1, 0))
## X-squared = 8.3768, df = 1, p-value = 0.0038

Model 29 — Model 1 but with no choice kernel

m29_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK21030.csv")
m29_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK89190.csv")
# Replace when back  
#m29_study1_eb <- m29_study1_eb_v1
m29_study1_eb <- rbind(m29_study1_eb_v1, m29_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m29_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK64103.csv")
m29_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK74167.csv")

# Delete when final one is back  
m29_study2_eb <- m29_study2_eb_v1
  
m29_study2_eb <- rbind(m29_study2_eb_v1, m29_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

# write.csv(m29_study1_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK_merged.csv")
# write.csv(m29_study2_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK_merged.csv")

Dramatically worse fit

ComparePars(m29_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m29", "m1")

ComparePars(m29_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m29", "m1")

sum(m29_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 1815.325
sum(m29_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 2777.038

Model Comparison

AIC_diff_df <-
 rbind(
    data.frame("diff"=m12_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Bayes"),  
    data.frame("diff"=m12_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Bayes"),
    
    data.frame("diff"=m4_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decayto0"),  
    data.frame("diff"=m4_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decayto0"),  
  
    data.frame("diff"=m3_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay"),  
    data.frame("diff"=m3_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay"), 
    
    data.frame("diff"=m5_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_beta"),  
    data.frame("diff"=m5_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_beta"),
    
    data.frame("diff"=m13_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Bayes_diff_beta"),  
    data.frame("diff"=m13_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Bayes_diff_beta"),
    
    data.frame("diff"=m6_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_lr"),  
    data.frame("diff"=m6_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_lr"),
    
    data.frame("diff"=m14_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_bayes"),  
    data.frame("diff"=m14_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_bayes"),
    
    data.frame("diff"=m15_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_diff_decay_bayes"),  
    data.frame("diff"=m15_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_diff_decay_bayes"),

    data.frame("diff"=m1_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_phi"),  
    data.frame("diff"=m1_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_phi")
 )

AIC_diff_df$study <- as.factor(AIC_diff_df$study)
AIC_diff_summ <- AIC_diff_df %>% group_by(model, study) %>% summarize(m=mean(diff))
## `summarise()` has grouped output by 'model'. You can override using the
## `.groups` argument.
sum(m1_study1_eb$AIC-m3_study1_eb$AIC)
## [1] -134.026
sum(m1_study2_eb$AIC-m3_study2_eb$AIC)
## [1] -147.2234

Confirm min in both studies is diff phi

AIC_diff_summ[AIC_diff_summ$study == 1, ] %>% arrange(m)
## # A tibble: 9 × 3
## # Groups:   model [9]
##   model              study     m
##   <chr>              <fct> <dbl>
## 1 Q_decay_diff_phi   1     -8.42
## 2 Q_decay_diff_lr    1     -7.36
## 3 Q_decay            1     -7.35
## 4 Q_diff_decay_bayes 1     -7.06
## 5 Q_decayto0         1     -6.38
## 6 Q_decay_diff_bayes 1     -4.83
## 7 Q_decay_diff_beta  1     -3.04
## 8 Bayes              1      4.98
## 9 Bayes_diff_beta    1      8.23
AIC_diff_summ[AIC_diff_summ$study == 2, ] %>% arrange(m)
## # A tibble: 9 × 3
## # Groups:   model [9]
##   model              study     m
##   <chr>              <fct> <dbl>
## 1 Q_decay_diff_phi   2     -9.95
## 2 Q_diff_decay_bayes 2     -9.03
## 3 Q_decay            2     -8.88
## 4 Q_decay_diff_lr    2     -8.69
## 5 Q_decay_diff_bayes 2     -6.42
## 6 Q_decayto0         2     -5.65
## 7 Q_decay_diff_beta  2     -3.74
## 8 Bayes              2      7.34
## 9 Bayes_diff_beta    2     11.8
levels <-
  c(
    "Bayes_diff_beta",
    "Bayes", 
    "Q_decayto0", 
    "Q_decay", 
    "Q_decay_diff_beta",  
    "Q_decay_diff_lr", 
    "Q_decay_diff_bayes", 
    "Q_diff_decay_bayes", 
    "Q_decay_diff_phi"
  )

assert(length(levels)==9)

AIC_diff_summ$model <- factor(AIC_diff_summ$model, levels = levels)
AIC_diff_df$model <- factor(AIC_diff_df$model, levels = levels)
mp1 <- ggplot(AIC_diff_summ %>% filter(study == 1), aes(x=model, y=m, group=model, fill=model)) + 
  geom_hline(yintercept = 0, size=2) +
  geom_jitter(data=AIC_diff_df %>% filter(study == 1), aes(x=model, y=diff, group=model), 
                fill="gray57", pch=21, width=.3, size=4, alpha=.05) + 
  geom_bar(stat="identity", color="black", alpha=.8, size=2) + 
  ylab(TeX('$\\Delta\ AIC$ from simple Q-learner')) +
  tol + ga + ap + tp +
  ggtitle("Study 1") +
  #facet_wrap(~ study) +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())  +
    theme(axis.title.y=element_text(size=22)) + ylim(-65, 65) + coord_flip()
mp2 <- ggplot(AIC_diff_summ %>% filter(study == 2), aes(x=model, y=m, group=model, fill=model)) + 
  geom_hline(yintercept = 0, size=2) +
  geom_jitter(data=AIC_diff_df %>% filter(study == 2), aes(x=model, y=diff, group=model), 
                fill="gray57", pch=21, width=.3, size=4, alpha=.05) + 
  geom_bar(stat="identity", color="black", alpha=.8, size=2) + 
  ylab(TeX('$\\Delta\ AIC$ from simple Q-learner')) +
  tol + ga + ap + tp +
  ggtitle("Study 2") +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())  +
  theme(axis.title.y=element_text(size=22)) + ylim(-65, 65) + coord_flip()
combined_aic <- mp1 + mp2
combined_aic

#ggsave("../paper/figs/fig_parts/aic_plot.png", combined_aic, height=5, width=10, dpi=700)
AIC_diff_df %>% group_by(model, study) %>% summarize(m=mean(diff)) %>% arrange(m)
## `summarise()` has grouped output by 'model'. You can override using the
## `.groups` argument.
## # A tibble: 18 × 3
## # Groups:   model [9]
##    model              study     m
##    <fct>              <fct> <dbl>
##  1 Q_decay_diff_phi   2     -9.95
##  2 Q_diff_decay_bayes 2     -9.03
##  3 Q_decay            2     -8.88
##  4 Q_decay_diff_lr    2     -8.69
##  5 Q_decay_diff_phi   1     -8.42
##  6 Q_decay_diff_lr    1     -7.36
##  7 Q_decay            1     -7.35
##  8 Q_diff_decay_bayes 1     -7.06
##  9 Q_decay_diff_bayes 2     -6.42
## 10 Q_decayto0         1     -6.38
## 11 Q_decayto0         2     -5.65
## 12 Q_decay_diff_bayes 1     -4.83
## 13 Q_decay_diff_beta  2     -3.74
## 14 Q_decay_diff_beta  1     -3.04
## 15 Bayes              1      4.98
## 16 Bayes              2      7.34
## 17 Bayes_diff_beta    1      8.23
## 18 Bayes_diff_beta    2     11.8

Get AICs for supplemental table

cat(round(sum(m1_study1_eb$AIC), 0), "\n")
## 42125
cat(round(sum(m1_study2_eb$AIC), 0))
## 42798
cat(round(sum(m3_study1_eb$AIC), 0), "\n")
## 42259
cat(round(sum(m3_study2_eb$AIC), 0))
## 42945
cat(round(sum(m4_study1_eb$AIC), 0), "\n")
## 42380
cat(round(sum(m4_study2_eb$AIC), 0))
## 43391
cat(round(sum(m5_study1_eb$AIC), 0), "\n")
## 42797
cat(round(sum(m5_study2_eb$AIC), 0))
## 43655
cat(round(sum(m6_study1_eb$AIC), 0), "\n")
## 42257
cat(round(sum(m6_study2_eb$AIC), 0))
## 42972
cat(round(sum(m13_study1_eb$AIC), 0), "\n")
## 44206
cat(round(sum(m13_study2_eb$AIC), 0))
## 45795
cat(round(sum(m14_study1_eb$AIC), 0), "\n")
## 42573
cat(round(sum(m14_study2_eb$AIC), 0))
## 43284
cat(round(sum(m15_study1_eb$AIC), 0), "\n")
## 42294
cat(round(sum(m15_study2_eb$AIC), 0))
## 42925
cat(round(sum(m11_study1_eb$AIC), 0), "\n")
## 43177
cat(round(sum(m11_study2_eb$AIC), 0))
## 44171
cat(round(sum(m12_study1_eb$AIC), 0), "\n")
## 43800
cat(round(sum(m12_study2_eb$AIC), 0))
## 45184